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ABSTRACT 

The problem of turbulent flow past a backward facing step is important 
in many technological applications and has been used as a standard test 
case to evaluate the performance of turbulence models in the prediction of 
separated flows. It is well known that the commonly used K-e (and K-0 
models of turbulence yield inaccurate predictions for the reattachment 
point in this problem. By an analysis of the mean vorticity transport 
equation, it will be argued that the intrinsically inaccurate prediction of 
normal Reynolds stress differences by the K-e and K-S models is a major 
contributor to this problem. Computations using a new nonlinear K-e model 
(which alleviates this deficiency) are made with the TEACH program. 
Comparisons are made between the improved results predicted by this 
nonlinear K-e model and those obtained from the linear K-e model as well as 
from second-order closure models. 
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residence at the Institute for Computer Applications in Science and Engi- 
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INTRODUCTION 


Among the various turbulence models in existence, the linear K-e 
and K-e, models are the most widely used by scientists and engineers to 
solve practical problems. The primary advantages of the models 
include their broad invariance properties and the relatively simple 
manner in which they can be incorporated into most existing 
Navier-Stokes computer codes which allow for a variable viscosity (see 
Speziale 1987). Furthermore, for unseparated turbulent boundary 
layers, these linear K-e and K-e models have been shown to provide 
excellent descriptions of the flow (see Rodi 1982). 

However, one major drawback of the linear K-e (and K-e) model is 
that it yields highly inaccurate predictions for the normal Reynolds 
stress differences. Consequently, for the types of turbulent flow 
where the normal Reynolds stress differences play an important role 
(such as secondary flows in a non-circular duct or separated flows), 
the linear models can give rise to considerable inaccuracies. There 
have been many efforts over the years to include nonlinear effects in 
the modeling of the Reynolds stresses within a two-equation format 
(see Lumley 1970, Launder and Ying 1971, Gessner and Emery 1976, and 
Saffman 1977). However, these models do not exhibit the general 
invariance necessary for the broadest range of application. Recently, 
a nonlinear K-e (and K-e.) model was obtained by making an asymptotic 
expansion subject to the constraints of dimensional and tensorial 
invariance, realizability, and material frame-indifference in the 
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limit of two-dimensional turbulence (see Speziale 1987). This model 
was shown to yield improved predictions for the normal Reynolds 
stresses in internal flows which are unseparated. However, for any 
proposed model to be sufficiently applicable, it has to be successful 
in predicting turbulent flows which have separated regions since these 
occur in a wide variety of problems which are of technological 
importance. 

In this study, the new nonlinear K-e model is incorporated into 
the TEACH computer code in order to analyze the problem of a 
two-dimensional turbulent flow past a backward facing step. Turbulent 
flow past a backward facing step has served as a primary benchmark for 
the performance of turbulence models in the prediction of separated 
flows (c.f., Abbott and Kline 1962; Briggs, Mellor and Yamada 1977; 
Kim, Kline and Johnston 1980; Eaton and Johnston 1981; Sindir 1982; 
Celenligil and Mellor 1985; and Chen 1985). The nonlinear K-e model 
will be shown to yield improved predictions for the reattachment point 
and turbulence intensities. Comparisons between the computed 
turbulence statistics and existing experimental data (along with the 
predictions of other turbulence models) will be made. In addition, 
the sensitivity of the results to the new empirical constant in the 
nonlinear model will be examined along with the effect of the Oldroyd 
derivative terms. The numerical results obtained will be discussed in 
detail along with other prospective future research. 
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2. THE PHYSICAL PROBLEM 

The problem to be considered in this study consists of the 
turbulent flow of an incompressible viscous fluid past a backward 
facing step (see Figure 1). The governing equations to be solved are 
the Reynolds equation and the continuity equation which are of the 
general form (c.f., Hinze 1975) 


, 3v 

p ( + rxx 
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VP+pVV+V*T 


where 


V-v = 0 

v 5 Mean velocity field 
P = Mean pressure field 
p = Density of the fluid 
p = Dynamic viscosity of the fluid 


( 1 ) 

( 2 ) 


and x is the Reynolds stress tensor whose components are given by 


where u is the fluctuating part of the velocity field. 

In order to achieve closure, equations relating the Reynolds 
stress tensor to the global history of the mean velocity field are 
needed. For the linear K-e model of turbulence, the Reynolds stress 
tensor is assumed to be of the form (see Hanjalic and Launder 1972) 
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are, respectively, the turbulent kinetic energy per unit mass, and the 

mean rate of strain tensor. Here, e is the dissipation rate of 

turbulence, and C is a dimensionless constant which assumes an 

V 

approximate value of 0.09. 

Typically, at high Reynolds numbers, the turbulent kinetic energy 
and dissipation rate are modeled by transport equations (c.f., 
Hanjalic and Launder 1972) of the form 
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where C^, C^, C^, and are dimensionless constants which take the 
approximate values of 0.11, 0.15, 1.43, and 1.92, respectively. 

The above set of equations are closed, and form the basis for the 
linear K-e model. Because of its simple structure, the K-e model can 
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be easily incorporated into any Navier-Stokes computer code which 
allows for a variable viscosity. This feature, together with its 
accurate predictions of thin turbulent shear flows (i.e., unseparated 
turbulent boundary layers), has made the model very popular with 
engineers and scientists (see Rodi 1982). Nevertheless, despite these 
advantages, the K-e model is known to yield highly inaccurate 
predictions for the normal Reynolds stresses. For example, in a fully 
developed turbulent channel flow (see Figure 2), the linear K-e model 
predicts that the normal Reynolds stresses are all equal, i.e., that 


T = X = T 

XX yy z z 


(9) 


which is in substantial contradiction of experimental data. To be 
more specific, the experimental data of Laufer (1951) for turbulent 
channel flow at a Reynolds number of 30,800 indicates that contrary to 
(9), 



where ||*|| denotes the maximum norm. 

Such serious errors in the normal Reynolds stress difference (x 
- x ) can yield significant inaccuracies in the calculation of a 

X X 

two-dimensional recirculating flow with the mean velocity 


v = u(x,y)i + v(x,y)i (11) 

A problem arises since such a flow is a solution of the vorticity 
transport equation (which determines the mean velocity) 
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where v = y/p is the kinematic viscosity and 


O) = 


9v 9u 


9x 9y 

is the mean vorticity. In the recirculation zone of turbulent flow 
past a backward facing step (see Figure 1), the term 


9 ( x - x ) 
v yy xx 7 
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(13) 


(which vanishes in a unidirectional mean turbulent flow or in an 
unseparated turbulent boundary layer) is of a comparable order of 
magnitude to the Reynolds shear stress terms on the right-hand-side of 
(12). This arises from the fact that x-deri vati ves are of a 
comparable order of magnitude to the y-deri vati ves and, as indicated 


above, i | t - t 
* 1 1 yy xx 


> || t | | in the recirculating zone. Hence, the 
xy 


serious inaccuracies that arise from the linear K-e model in the 
calculation of the reattachment point for the backward facing step 
problem may be largely due to this modeling deficiency in the normal 
Reynolds stresses (see Speziale 1987). 

The linear K-e model substantially underpredicts the reattachment 
point for turbulent flow past a backward facing step as discussed 
extensively at the 1980-81 AFOSR-HTTM Stanford Conference on 
Turbulence (the separation length L/aH predicted by this model varies 
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from 5.2 to 5.5, whereas experiments indicate that it should be 
approximately 7.0). As alluded to above, a significant cause of this 
problem may arise from the inability of the linear K-e model to 
predict normal Reynolds stress differences accurately. Recently, 
Speziale (1987) derived a nonlinear generalization of the K-e model 
which takes the form 
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and Cp = C^. = 1.68. As a result of the presence of the Oldroyd 

o _ 

derivative term D . j and the quadratic terms in D. j , this nonlinear K-e 
model is able to describe turbulent memory effects (it bears a certain 
resemblance to the Ri vl in-Ericksen fluids of viscoelastic 
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where (for the nonlinear K-z model) 
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In deriving (18) - (21), the continuity equation (2) given by 
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( 22 ) 


has been made use of. Equations (16) - (22) must be supplemented with 
transport equations for K and £. For the problem under consideration, 
turbulent diffusive effects can be neglected (see Chen 1985) and, 
thus, the transport equations (7) - (8) can be simplified to the form 
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Equations (16) - (24) represent a closed system of equations for the 
determination of the mean turbulent flow fields. Of course, these 
equations must be solved subject to the appropriate boundary 
conditions. For turbulent flow past a backward facing step, these 
boundary conditions consist of fully-developed turbulent channel flow 
sufficiently far upstream and downstream of the step. The law of the 
wall is applied at the solid boundaries (c.f., Rodi 1982). More 
details concerning the boundary conditions will be presented in the 
next section where the numerical approach will be discussed. 



11 


3. NUMERICAL METHOD 

In this study, the Imperial College's TEACH (Teaching Elliptic 
Axisymmetry Characteristics Heuri stically) computer program will be 
used since it was written for the calculation of recirculating flows 
and performs reasonably well for such problems (see Lilley and Rhode 
1982). Furthermore, because this program is based on the linear K-e 
model, it is a relatively straightforward process to incorporate the 
nonlinear K-e model (this simply requires the addition of the 
appropriate nonlinear terms where the Reynolds stresses appear). 

The TEACH computer program will be implemented on the nonuniform 
finite difference mesh shown in Fig. 3. This 166x73 finite difference 
mesh spans 5 step heights upstream and 30 step heights downstream from 
the step. The grids are uniform and extremely dense in the 
recirculation region and in the vicinity of step corners in order to 
conveniently calculate and adequately resolve the higher-order 
velocity gradients. Special merits of the staggered grid system used 
in the TEACH code are discussed in Gosman and Pun 1974, and Gosman and 
Ideriah 1976. The scheme is second-order accurate overall since the 
first-order derivatives are evaluated by central differences over a 
single mesh spacing (c.f., Roache 1972). However, in evaluating 
second-order derivatives (which are needed in the nonlinear model), 
this nonuniform mesh requires a more complex finite difference 
expression which has a second-order accurate form given by 
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where f(i) is any function, f 1 ( i ) is its first derivative, f"(i) is 
its second derivative, A + is the mesh spacing in the forward 
direction, and A_ is the mesh spacing in the backward direction. 
Equation (25) can be applied for the calculation of second partial 
derivatives with respect to x or y (by invoking the continuity 
equation, no mixed derivatives in x and y need to be calculated). 
Since this finite-difference expression requires the value for f l , 
which is not exact, this expression is actually somewhat less than 
second-order accurate. Hence, the application of equation (25) is 
limited to outside of the recirculation zone where second derivative 
contributions are small. Within, and immediately adjacent to, the 
recirculation zone (where the mesh is uniform), fourth-order accurate 
expressions for the first and second derivative terms in the Reynolds 
stresses are used that are given by 
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A top hat filtering scheme is used to smooth any numerical 

fluctuations which arise due to the— steep velocity gradients. The 

reader is referred to Isaacson and Keller (1966) for more details on 
these numerical representations. 

In the TEACH code, the governing equations for the linear K-e 
model are the continuity equation, the x and y components of the 
Reynolds equation, and the kinetic energy and dissipation rate 
transport equations as given in the previous section. The Poisson 
equation for pressure (obtained by taking the divergence of the 
Reynolds equation (1)) is also used since the solution procedure here 
is based on primitive variables (i.e., is based on pressure-vel oci ty 
as opposed to the stream function-vorti ci ty approach). Since the 
problem considered is at very high Reynolds numbers (above 100,000), 
the contributions of the diffusion terms in the transport equations 
for K and e are quite small and can be neglected. For more details on 
this general numerical method, the reader is referred to Lilley and 
Rhode (1982). 

The governing equations for the nonlinear K-e model are the same 
as for the linear model except for the differences in the expressions 
for the Reynolds stresses. Since the additional nonlinear terms in 
the Reynolds stresses are not extremely large, they are simply treated 
as added source terms in the TEACH code. It should be noted that the 


vari abl es 
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are already calculated for the linear K-e model and, hence, are 
generated by the basic TEACH code. The only new variables that need 
to be calculated are the second derivatives of the velocity fields 
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Inside the reci rcul ation region, these second derivatives are computed 
by the fourth-order accurate finite difference scheme given by 
equation (27). A lower-order accurate scheme was tried, but it 
yielded some fluctuations in these derivative fields (especially 
within the recirculation zone) due to its inability to resolve the 
steep velocity gradients. This problem persisted (although to a 
substantially lesser extent) even with the higher-order accurate 
scheme and therefore a smoothing filter was employed. As alluded to 
earlier, a top hat filter was used. This filter can be applied to a 
field more than one time, yielding a smoother result after each pass. 
The velocity fields are filtered once before they are used in 
evaluating first and second derivatives. As mentioned before, a 
fourth-order accurate expression is also used to calculate the first 
derivatives within the recirculation zone in the nonlinear Reynolds 
stress calculations. It should be noted that at the solid boundaries, 
an outward Taylor expansion is used to evaluate the derivatives. 
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After obtaining the first and second order derivatives, the nonlinear 
part of each Reynolds stress is readily obtained from equations (18) - 
(21). In order to incorporate these nonlinear Reynolds stress terms 
into the program, the linear representations for the Reynolds stresses 
were located in the various parts of the TEACH code, and the nonlinear 
terms were added in an appropriate fashion. 

Computations were conducted in a channel with an expansion ratio 
of 3:2 at a Reynolds number of approximately 132,000 (based on the 
upstream centerline mean velocity and downstream channel width). The 
program was run with an automatic vectorizer on the CYBER 205 computer 
at the University of Georgia. Iterations were performed until a 
converged solution was obtained based on a residual source criterion 
(see Lilley and Rhode 1982). Approximately 90 minutes of CPU time 
were required on a 166x73 mesh for the nonlinear K-e model to 
converge, whereas only about 45 minutes were needed for the linear K-e 
model. The number of steps required for a converged solution depends 
on the model and the input conditions and varied from 400 to 1000 
iterations. The closer that the initial guess of the variable fields 
is to the actual solution, the less iterations are needed. However, 
by a comparison of the rate that the residual sources decrease, it was 
clear that the nonlinear model converges at a slower rate than the 
linear model. This is not very surprising considering the fact that 
the nonlinear K-e model contains Reynolds stress relaxation terms 
which are dispersive rather than dissipative. The detailed numerical 
results obtained will be discussed in the next section. 
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4. DISCUSSION OF THE RESULTS 

We will now present the computed results obtained for turbulent 
flow over a backward facing step for a Reynolds number Re = 132,000 
and an expansion ratio of 3:2. The computed streamlines and mean 
velocity profiles obtained from the linear K-e model are shown in Fig. 
4 which clearly indicate a reattachment length of L/aH = 5.5 -- a 
value in the range of previously conducted computations. This 
compares rather unfavorably to the experimental value of L/aH = 7.0 
(see Kim, Kline, and Johnston 1980). The nondimensional turbulence 
intensity (uu) and shear stress uv obtained from the linear K-e 
model are shown in Figs. 5-6 alongside the available experimental data 
(reliable experimental data is not available for the recirculation 
zone). The computed streamlines and mean velocity profiles obtained 
from the nonlinear K-e model, shown in Fig. 7, clearly demonstrate an 
improved prediction for the reattachment point of L/AH = 6.4. This 
improvement is probably due to the better prediction of the turbulence 
intensities in the recirculation zone (the reader should compare Figs. 
8-9 with Figs. 5-6). However, since no reliable Reynolds stress data 
is available inside the recirculation zone (due to flow oscillations), 
the present conclusion must rest on comparisons of data from Kim, 
Kline, and Johnston (1980) downstream of the reattachment point. 

Comparisons will now be made between the results of the nonlinear 
K-e model and those obtained from second-order closure models. As 
shown in Fig. 10, the second-order closure model of Briggs, Mellor, 



17 


and Yamada (1977) predicts a reattachment length of L/AH = 8.0 (it 
should be noted, however, that in a more recent study of this model 
conducted by Celenligil and Mellor 1985, a result L/aH = 7.7 was 
claimed which is of a comparable accuracy to the result of L/aH = 6.4 
predicted by the nonlinear K-e model). The reasons for this 
compatibility are more apparent when Reynolds stress results are 
compared. As shown in Fig. 11, the turbulence intensity predictions 
of the nonlinear K-e model and the second-order closure model are in 
good qualitative agreement throughout the unseparated flow. The 
turbulent shear stress predictions shown in Fig. 12 are not in as good 
agreement inside the recirculation zone (it is unfortunate that no 
reliable data is available in this region in order to make a critical 
comparison). However, in so far as the reattachment point is 
concerned, it is clear that the nonlinear K-e model yields comparably 
good, if not better, results to second-order closure models with the 
need for substantially less calculations (second-order closure models 
require a higher level of computation since transport equations must 
be solved for each component of the Reynolds stress tensor). 

The nonlinear K-e model examined herein introduces two new 
empirical constants, Cq and C^, in equation (14). These constants 
were found to assume a value of 1.68 by correlating with experimental 
data for normal Reynolds stress differences in turbulent channel flow 
(Speziale 1986). However, the accuracy of this data is somewhat 
questionable; errors of the order of 10% can easily occur. 
Consequently, calculations were performed to test the sensitivity of 
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the computed results to the precise value of Cp and C^. For this 
purpose, the computations were repeated for Cp = = 1.40. The 
computed streamlines and turbulence intensity for this case are 
compared with the previous computations in Figs. 13 - 14. There is no 
significant difference in the two results, indicating that even a 15% 
error in the predicted value of Cp and would have little effect on 
the major conclusions of this study. 

Computer runs were also made for the nonlinear K-e model where 
the second-order derivative terms of the mean velocity arising from 
the Oldroyd derivative in (14) are set to zero (such a resulting model 
bears a qualitative resemblance to previous nonlinear models proposed 
by Lumley 1970 and Saffman 1977). As shown in Fig. 15, without the 
second-order derivatives, the resulting reattachment length reduces a 
significant amount to a value of L/aH = 6.0. This shortening of L/aH 
is not unexpected since the second derivative terms have a dispersive 
character which would reduce the dissipation in the separated zone 
allowing it to expand. It is rather surprising, however, that the 
elimination of the second derivative terms had a relatively small 
effect on the turbulent stress intensity shown in Fig. 16. Since the 
second derivative terms generally are of a comparable order of 
magnitude to the first derivative terms squared (see equations (18) - 
(21)), one might expect the contribution of these terms to be more 
significant in the turbulence intensity predictions. Hence 
considering their substantial impact on the reattachment point, the 
most dominant contribution of the second derivative terms must be in 
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the prediction of normal Reynolds stress differences (quantities which 
play an important role in the calculation of a recirculating flow as 
discussed earlier). Thus, it is clear that the Oldroyd derivative 
(which distinguishes this new nonlinear K-e model from all previous 
such nonlinear models) plays a significant role in the turbulent 
backward facing step calculation. 

Unfortunately, as discussed in Section 3, the accurate evaluation 
of higher-order derivatives is extremely difficult even with a 
fourth-order accurate finite difference scheme. In this study, a top 
hat filter was used to alleviate some localized numerical fluctuations 
in these derivatives. Since the top hat filter is an averaging 
algorithm, it also somewhat reduces the steep velocity gradients in 
the flow, thus adding to the problem of numerical accuracy. However, 
this problem is a localized one and any adverse consequences of 
filtering appears to be in the direction of underpredicting the 
reattachment point as a result of smoothing. Some initial 
computations that were conducted without filtering (which must be 
viewed with suspicion as a result of fluctuations in the second 
derivative terms) yielded a reattachment length of L/aH = 6.7. Thus, 
it appears that the more accurate calculation of the second 
derivatives is likely to bring the computed results of this study in 
closer agreement with experimental observations. It is our opinion, 
however, that such improvements in the numerical algorithm are 
unlikely to make more than a 5 or 6% change in the numerical results 
presented herein. 



5. CONCLUSIONS 


The nonlinear K-e model examined in this study has been shown to 
yield considerably improved predictions for the reattachment point for 
turbulent flow past a backward facing step. To be specific, the 
nonlinear K-e model was shown to predict a separation length of 
L/aH = 6.4 as compared to the experimental value of L/aH = 7.0. This 
result constitutes a substantial improvement on the value of L/aH in 
the range of 5.2 - 5.5 predicted by the linear K-e model and is of a 
comparable accuracy to results obtained from the substantially more 
complicated second-order closure models (e.g., Celenligil and Mellor 
1985 obtained a separation length of L/aH = 7.7). From an analysis of 
the mean vorticity transport equation, it was argued that these 
improved results probably arise from the nonlinear K-e model’s ability 
to predict normal Reynolds stress differences more accurately. While 
there is no reliable experimental data for the Reynolds stresses 
inside the recirculation zone, the results obtained from the nonlinear 
K-e model are in the range of those obtained from second-order closure 
models and appear to constitute an improvement over those predicted by 
the linear K-e model . 

Several questions about the accuracy of the numerical results 
presented in this study still need to be resolved as is true of most 
numerical studies of the backward facing step problem. The TEACH code 
has certain undesirable features in its convergence properties when 
terms (such as the nonlinear contributions to the Reynolds stresses) 
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which are not purely dissipative are added. Furthermore, as a result 
of the large velocity gradients in the recirculation region, it was 
extremely difficult to calculate higher order derivatives of the 
velocity field in a highly accurate fashion. However, as a direct 
consequence of the smoothing properties of the top hat filter, it 
appears that any errors in the computed reattachment point would be in 
the direction of underprediction, thus, putting the results obtained 
herein in closer agreement with the experimental data. Nonetheless, 
it would be useful to check these results using several alternative 
numerical algorithms. Likewise, it would be of value to consider 
other geometries. Unfortunately, such investigations are beyond the 
scope of the present study and must await future research. Although 
several questions remain to be answered, the results obtained in this 
study are extremely encouraging and strongly support the pursuit of 
future investigations of this nonlinear K-e model which may prove to 
be of considerable value in the future analysis of a variety of 
turbulence problems. 
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FIGURE 2. Fully— developed turbulent channel flow. 
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Dimensionless turbulence intensities obtained from the linear K — e model. 
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Turbulence intensities obtained from the nonlinear K — e model 
(a) C D = C E = 1.40 and (b) C 0 = C E = 1.68 
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FIGURE 15. Computed streamlines obtained for the nonlinear K-e model 
without the second derivative terms. 
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